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Recovery of the sound speed for the Acoustic 
wave equation from phaseless measurements 


Joonas Ilmavirta* Alden Waters^ 


Abstract 

We recover the higher order terms for the acoustic wave equation from 
measurements of the modulus of the solution. The recovery of these coeffi¬ 
cients is reduced to a question of stability for inverting a Hamiltonian flow 
transform, not the geodesic X-ray transform encountered in other inverse 
boundary problems like the determination of conformal factors. We obtain 
new stability results for the Hamiltonian flow transform, which allow re¬ 
covery of the higher order terms. Previous techniques do not measure the 
full amplitude of the outgoing scattered wave, which is the main novelty 
in our approach. 


1 Introduction 

Scattering is a general physical process where some forms of radiation, such as 
light or sound, or moving particles are forced to deviate from a geodesic trajec¬ 
tory by a path due to localized non-uniformities in the medium through which 
they pass. Conventionally, this also includes the deviation of reflected radiation 
from the angle predicted by the law of reflection m- Scattering may also refer 
to particle-particle collisions between molecules, atoms, electrons, photons and 
other particles m- The types of non-uniformities which can cause scattering 
are sometimes known as scatterers. These include particles and surface rough¬ 
ness. The type of stability estimate we prove on the phaseless measurements 
of the solution insure that the energy of the waves uniquely determines the 
scatterer. 

The theory of signal processing and inverse problems has seen a recent in¬ 
crease in a class of so-called phaseless measurements. Often in experiments, 
when a source wave is measured, the only part of the information available to 
experimenters is the modulus of the wave from the source. In signal processing, 
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algorithms found in [T^ and m are focused on the recovery of waves from a 
sequence of Fourier modes. We are interested in phaseless measurements to re¬ 
cover scattering terms for the acoustic wave equation and generalized Helmholtz 
equation. The modulus of the wave corresponds to the energy density of the 
wave at a given point. The inverse boundary value problem differs from the 
questions examined in signal processing where often times one is dealing with 
incomplete data sets. In particular, we show that an idea of which partial dif¬ 
ferential equation the wave comes from is enough to give a full reconstruction of 
the coefficients modulo diffeomorphism. The scattered wave then completely de¬ 
termines the scatterer. These results are supported by the numerical work in [5] 
and are applicable to other operators which admit a Gaussian beam type solu¬ 
tion. The problem differs from the both of the author’s previous work inn [ID 
because the terms which we are recovering come from higher order terms which 
control the bicharacteristic flow associated to the Hamiltonian governing the 
partial differential equation. 

Practical applications to phaseless problems are varied — and one such ap¬ 
plication is multi-wave tomography. In multi-wave tomography usually some 
type of wave is sent to a portion of the body which is being imaged. In electro¬ 
magnetic or optical radiation tomography the wave interaction with the tissues 
of the patient are measured (D. Naturally one cannot measure inside the pa¬ 
tient, so some initial boundary value problem must be considered. Similarly, 
to image the Earth, one has to send waves of some kind through the planet 
and make measurements at the surface. One such mathematical model of the 
emitted ultrasound waves is the acoustic wave equation with a high-frequency 
source term. 

Let M C be a bounded and smooth manifold. Let g be a Riemannian 
metric on which agrees with the Euclidean one outside M. Let the standard 
Laplace-Beltrami operator be denoted as 

in local coordinates with g{x) = {gik{x)), and (g®*(x)) = {gki{x))~^. We con¬ 
sider manifolds, M, which are smooth {C°°). We write the local coordinates as 
(a;^, ..,x'^). We also assume the manifolds have a boundary. 

The generalized Helmholtz equation may be written as 

Lu = AgU + {iaX + X^)n^{x)u = h{x, X) xGM.‘^. (1.2) 

The scalar A is large and n?(x) is sound speed. The number a helps describe 
the attenuation coefficient and a € (0,1). It is zero in the limit of zero absorp¬ 
tion. The source h(x, A), which emits the waves, we also assume to depend on A 
and be compactly supported in x with codimension 1. Sources will be modeled 
after Gaussian beams following [32) . We may pick sources anywhere inside the 
domain M, however we chose a particular set of them in order to provide a com¬ 
plete reconstruction of n^(x). While this equation is known as the generalized 
Helmholtz equation, multiplication by a prefactor, exp(*At), takes the solution 
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u and turns it into a solution exp{iXt)u to the acoustic wave equation when 
n^{x) = 1 on dM, and a —?► 0 

dfU — n~^{x)AgU = exp{iXt)h{x, X) x S (1.3) 

The acoustic wave equation models the scattering of waves in the Earth’s core. 

We chose to model our solutions to the wave equation with Gaussian beams. 
The existence of Gaussian beam solutions to the wave equation has been known 
since the 1960’s first in connection with lasers, see Babic and Buldyrev [7]. 
They were also used in the analysis of propagation of singularities in PDEs by 
Hormander m and Ralston [?D]. In the context of the Schrodinger equation 
first order beams correspond to the so-called classical coherent states. The 
higher order versions of these equations have been introduced to approximate the 
Schrodinger equation solutions in quantum chemistry by Heller [18], Hagedorn 
and Herman and Kluk [19]. 

Given a smooth, strictly convex, bounded domain M equipped with a metric 
g we assume that on the boundary of M, dM that n^(x) = 1. The measurements 
we consider give then data of the form 


{(x, |u(x)|) : X e dM} 


(1.4) 


with the collection of {x, X) varying over all Xq G dM and ujq G S^M, 

which are inward pointing into the manifold. We need this collection of data 
in order to give a complete reconstruction of n^(x). We consider the metric 
g to be fixed and n^(x) to vary. This collection of measurements is a ‘true’ 
phaseless problem in contrast to phaseless backscattering measurements which 
were recently investigated in [27]. We see that the measurements for both 
equations (1.31 and (1.2) coincide. Thus we have developed a robust model 


for the scattering of both waves traveling through the earth’s core (1.3) and 


quantum mechanical particles in (1.2). 


In this paper, we derive a stability result for the higher order coefficients 
of the acoustic wave equation (1.2) for fully phaseless measurements. Stable 


reconstructions have been made from Robin conditions for lower order terms 
than considered here [3]. In the related case, for the generalized Helmholtz 
equation for Dirichlet boundary conditions in [241 135] and Robin conditions 
in [a EH] the potential can also be recovered. However stability estimates from 
phaseless measurements have not been previously given. In [29[ [28] , uniqueness 
results in dimension 3 for lower order terms than the ones considered are de¬ 
rived from phaseless measurements. These papers are predicated on analyticity 
arguments, which require data in a small neighborhood of the source. The ques¬ 
tion of phaseless stability from internal measurements for Schrodinger was also 
examined in [1] . These measurements are in contrast to the boundary data we 
require. We are not able to prove uniqueness results unless A —>■ oo. However, 
physically this is equivalent to setting the n^(x) term equal to zero. The inverse 
problem of recovering the source for wave equations from Dirichlet boundary 
conditions is also examined in [3H1 Ell HOI ESI Ei] ■ The major difference is that 
we are able to recover higher order terms which make the bicharacterstic flow 
different. 


3 







2 Statement of the Main Theorem 

Our main theorem can be stated as follows. We consider sources X). 

with support in Here we view M as an embedded submanifold of M', 

a larger manifold. Sources are indexed by xq and wq. We let v{x) denote the 
outward unit normal to the boundary at x. Let the set 

d+SM = {(xo,wo) : xq G dM, {iy(xo),u)o) > 0, wq G Sa:M}, (2.1) 

be the range of Xq, Wq where S^M denotes the unit sphere bundle of the manifold 
at X. The sources are defined for x G as: 

A) = 2iXx\{x) (wo • i' + Ad(0)(x — xq) • i' + VgXA • v) exp{iXtjj{x)) 

( 2 . 2 ) 


with 

v>(x) = (x - Xo) • Wo + ^(a; - Xo) • M{0){x - xq). (2.3) 

Here X\i^) is a smooth compactly supported function of codimension 1 in a 
neighborhood of the set 

{x G M', |x — xol < (2.4) 


which is perpendicular to wq. This set also contains the support of 'ipix). The 
flow for the ray path of H is defined by the set of ODEs in local coordinates on 
the manifold: 


^ - 0 _ ^g-^^(a^(s))PfcPi dn^{x{s)) 

ds ^ X s pj, 


(2.5) 


with initial data (xo,wo). The matrix A4(0) is a complex matrix satisfying 

A4(0) = A4(0)^, A4(0)x(0) = p(0), 9A4(0) positive definite on x(O)'''. 

( 2 . 6 ) 


We know there is a solution to the following equation for each xo,wo: 

(Ag + (A^ + iaA)n2(x))u"»’“» = ^“’‘"“(x, A), (2.7) 

We let eo G (0,1) and fi^ = n\ — n|. We have the following stability result 
on where the measurements are the amplitudes of the collection of 
corresponding to different n^, ranging over sources indexed by d+SM. Let 
diam//(M) denote the maximal radius of the manifold with respect to the flow 
H. We let g' denote the extended metric which is (/ on M and Euclidean on the 
exterior. We impose the condition the manifold M' contains M and be such 
that M' is simple with respect to the metric n~^g'. We recall that a simple 
manifold is one which is strictly geodesically convex with respect to the metric 
g, and has no conjugate points. 
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Theorem 1. Let N > max{(l — d)/2 + 2 + 2s, 0} + 1, s > d/2. Then there 
exists a constant Ci, which depends on diam//(M), the C'^{M) norm of , 
and a constant C 2 , which depends on diain//(M) and the norm of 

n'j{x),i = 1,2, such that if and solve the radiation problem with 

coefficients n\ and n\ then it follows that if 

|Pl - «2|U3fM^ < ^ < ^0’ °= IPl I - 1^2 II < eo, 

xGdM;{xo,ujQ)Gdj^SAI 

( 2 . 8 ) 

then this implies 

11’^2 ~ ''^1 1 Il2(m) — 

for some /3' € (0,1) depending on dimension d, with Ae > 1. 

The uniqueness corollary follows immediately. 

Corollary 1. Assume d = 0 and that the assumptions of Theorem^are satisfied 
for all large A. Then n\ = n\. 


Our results indicate that there is very little stability to be expected from 
such a problem. 

Remark 2. The standard problem see for example JdtA \ 34^ is to consider 

(-n-^Ag - A2 )u = 0, (2.10) 

u\aM = h{x,\). 


with h{x, A) a compactly supported function on dM. Here one would recover the 
factor from the Dirichlet-to-Neumman map with h{x, A) a high frequency 
source, with measurements of the form 

f \d,.iui - U2)f dgS (2.11) 

dM 


We see that 


[ \\u,\^-\U2\^\dgS 
dM 


< 


Ul 


U2\'^dgS 


( 2 . 12 ) 


dM 


with ui,U 2 corresponding to different nt. Here dgS denotes the surface measure 
on M. Our problem is more physical because it corresponds to measuring the 


amplitudes of two different waves as on the left hand side of (2.12), rather than 


super-imposing them and measuring their difference. We are also able to include 
an attenuation coefficient which is the a —> 0 limit of (2.10). 


5 




Generically, we are looking for an asymptotic model to (1.2) of the form 


U{x) = A{x, A) ejq){iX'ip{x)) = ^ ex.p{iX-ip{x)), (2-13) 


which we show in the high frequency limit solves the equation (1.2) up to suitable 
error terms. The variable x G M and s = s{x) is a parameter which helps 
describe the bicharacterstic flow in terms of the coordinates on M. We use a 
Gaussian beam Ansatz which involves the construction of a phase function tp^x) 
and an amplitude a for the zeroth order beams in local coordinates as 


tp{x) = S{s) + (a; — x{s)) ■ p{s) + — x{s)) ■ M{s){x — x(s)) 

a{s,x) = ao{s) + 0{dg'{x,x{s))) 


(2.14) 


where x(s) is a curve which describes a Hamiltonian flow, given by (2.5) below, 
and S'(s), Al(s) will be specified functions of s determined by the Gaussian beam 
Ansatz. For the construction of the Ansatz, we follow the work of [35] quite 
closely. In the previous work we considered the operator 


while here we use 


L — A]]^d T A iXoLTi?'{x') j 


L = Ag + X^n^(x) + iXav?[x). 


(2.15) 


(2.16) 


This corresponds to the acoustic wave equation. These operators have cor¬ 
responding Hamiltonian flows given by the Hamiltonian functions which are 
H = IpP — 1 and H — |p|g — n^. The main problem of isolating the X-ray 
transform of the coefficients n^(x) is to control the null-bicharacteristics. In 
local coordinates, the ordinary differential equations which govern the ray path 
{(s,x(s)) : 0 < s < T} along which solutions are concentrated are given by (2.51 
which is associated to L, resp. 


dPx{s) 

ds"^ 


(2.17) 


which is associated to L. The second equation gives that the ray paths in 
are straight lines while the first one does not. The problem of recovering n^(x) 
for the operator L is more difficult, and we address it in this paper. 


In the Galderon Problem in Gonformally Transversal Geometries, [26] and 
also [H] reduces the question of boundary distance rigidity to a question of 
invertibility of the geodesic X-ray transform. As a consequence of their work, 
they reduce the question of recovery of source terms for several operators from 
the Dirichlet-to-Neumann maps to a question of invertibility of the geodesic 
X-ray transform. In this paper, we choose phaseless data as our measurements 
and show that the question of recovery of sound speed amounts to a question 
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of invertibility and stability of a so-called flow transform. The results presented 
here have applications to other operators for which the Hamiltonian flow and 
the geodesic flow do not coincide. In order to prove our results, we intro¬ 
duce a condition which one can think of as a generalization of the condition of 
Bardos-Lebeau-Rauch [30] for Hamiltonian flows. Specifically, we require the 
Hamiltonian flow to be simple. 

Notation: For two functions f,g, we write / ^ g if there exists a constant 
C > 0 such that C~^f < g < Cg. We denote dg{x,y) the distance between the 
points cc, y defined by the Riemannian metric g. 


3 Construction of Solutions 

We make more precise the explicit nature of the solutions. We will use the 
Ansatz, 


i 

U{x) = exp{iXtf:{x))aj{x)X~^ 
j=o 


(3.1) 


to build asymptotic solutions to (1.2) in the high frequency limit, and then give 
estimates on the difference between these approximate solutions and the true 
solutions. We will follow [33] for the Gaussian beam Ansatz in constructing the 
phase a{x) and amplitude We claim: 


Theorem 3. There is a zeroth order Gaussian beam which solves the Helmholtz 
equation (1.2) in the free space with no source term. In local coordinates it takes 
the form 


U{x) = {ao{s) + 0{dg{x,x{s)))ey:p{i\il:{x))4i{x) + 0{X ^), (3.2) 

where x{s) is a curve in space defined by ( |2.5| ) with initial condition (xo,ujo) € 
d+SA4 and (p{x) is a cutoff which will be defined in a neighborhood of the curve 


ao(s) 




(3.3) 


and ip{x) is given by (2.14) and A4{s) is a positive definite complex matrix 
whose evolution is governed by a Ricatti equation. 


The A^*^ order Gaussian beam Ansatz can be constructed accordingly. How¬ 
ever, for the purposes of this paper the zeroth order terms are the most impor¬ 
tant. 


Remark 4. One might be tempted to use a real-phase Ansatz which should be 
of the form 


exp(iA'0(a:)) 



\x -Xo\J ’ 


(3.4) 
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with s a parameter describing the bicharacteristic flow, and '4){x) a purely real 
function, A{x) encodes information about the flow transform of n^{x) or a sim¬ 
ilar solution found in \2,^ . However, it is not possible to obtain information 
about the coefficient without the use of a weighted Sobolev space and this is 
more difficult since this clearly introduces singularities in the denominator. The 
Gaussian beam Ansatz here was chosen because of the good control over the 
error estimates in higher Sobolev norms, which allows for the use of the embed¬ 
ding theorems in the section on observability inequalities. Also, because of the 
presence of the tail, one does not need to construct solutions which are localized 
along a central ray. This is necessary for the estimates (8.6) and (8.18) to hold. 


since Green’s theorem is no longer available. However, the real phase construc¬ 
tion has been efficient in recovery of coefficients for the Dirichlet-to-Neumann 
problem see 


We need the following lemma to start off with: 


Lemma 1. Every geometric optics solution concentrates on an open set around 
the ray path {(s,a:(s)) : 0 < s < T} for some finite T, with T > 0. The flow 
for the ray path of H is defined by the set of ODEs in local coordinates on the 
manifold: 


^ cIk dg^^{x{s))pkPj , dn^{x{s)) 

ds ^ ’ ds dxi dxi 

and also 

^ifixis)) = VgV'(a;(s)) ■ NpH{x{s),pis)) = n^(x(s)) 


(3.5) 


(3.6) 


Proof. We have H(x,p) = \p\^g — associated to H{x,Ngif{x)), so we set 
p = Ng'if{x). We are looking for a solution to H{x,Ng'tf{x)) = 0 because we 
want to solve ( |3.12| ) to high order. We differentiate H{x,Vgip{x)) with respect 
to X. This gives the following relationship: 


\/,H{x, Vgf:{x)) + D'^fj{x)VpH{x, Ngf:{x)) = 0 (3.7) 

where represents the Hessian associated to if. For any curve y{s) we have 
the following identity 


^Nifiyis)) = D\if{y{s))^ = (3.8) 

D^{iP{y{s)) (^-NpH{y{s),Ngif{y{s))^-N,H{y,Ngif{y{s)) 


by substitution of (3.7). It follows that if we are seeking a null-bicharacteristic 
curve that 


^ = NpHiixis)), Vgifixis))) (3.9) 

= ^(VgV'(a;(s)) = -N,,H{x{s),Ngif{x{s))) 

8 













so that (3.81 vanishes. Moreover we see for our particular H that 


ds 


i^{x{s)) = \/gilj{x{s)) ■ \/pH{x{s),p{s)) 


n^{x{s)). 


(3.10) 

□ 


Proof of Theorem We have that 

i 

LU = ^ exp{iXtf{x))cj{x)X~T (3.11) 

i =-2 

The coefficients Cj, j = 0,1,..., are defined recursively as follows 

c _2 = (n'^ix) - \\/gtp\'^)ao{x) = E{x)ao, 

c_i = iavf{x)ao + Vg • {ao\7gtf) + VgOo • Vg')/’ + E{x)ai, 

Cj = ian'^{x)ag + l + Vg • {aj + lVgfj) + VgOg + l • Vg'0 + E{x)aj + l + AgUj. 

(3.12) 

If we Taylor expand the coefficients a.j{x) around the central ray cc(s), we are 
looking for a system of differential equations for the coefficients. If we have 
chosen c _2 and c_i to vanish on the ray to third and first order respectively 
and define S as in |2.14[ This leads to the following set of differential equations: 

dsS{s) = 2nf{x{s)), dsao{s) =—tT{Ai{s))ao — an^{x{s))ao. (3.13) 

Denoting Mjk(s) = dl.^^'tp{x{s)), we also must have 

-dsMjkis) = {dl^p^H{x{s),p{s)))Mjt{s)Mkr{s) (3.14) 

+ H (^x,P,^( 2 ;(s),p(s))'^^fe(s) +5x^p,^(2;(s),p(s))M^g(s)) (3.15) 

I 

+dYkH{x{s),p{s)) (3.16) 

with H{x,p) = p • Vpi/ {x,p). 

These equations correspond to differentiating the equation 

dsi^ix)+XJg'ii){x) ■ H{x,XJgif) =Q (3.17) 

with respect to x and evaluating along the curve x{s) with respect to the con¬ 
straint El{x, ygif) = 0. This equation is the correct one to consider by the 
derivation of ( |3.6[ ). These equations can be easily checked |12]. The matrix 
A4 is known as the Hessian matrix. Higher order beams can be constructed 
by requiring c _2 vanishing to higher order on 7 . One can require the c's with 
j > —2 also vanish to higher order and obtain a recursive set of linear equa¬ 
tions for the partial derivatives of oq, ai,..., a;. More precisely, for an order 
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beam I = [iV/2] — 1, and Cj{x) should vanish to order N — 2j — 2 whenever 
—2 < j < I — 1. We construct higher order beams, but we are only inter¬ 
ested in the precise form of the zeroth order terms as they are sufficient for our 
asymptotics. We see that 


00(5) = X\{x) exp I J —an'^{x{t)) — trA^(t) 


dt 


The phase ip needs to verify the conditions 

ip{x{s)) = S{s), V^/>(x(s)) =p(s), D‘^iI;{x{s)) = M{s), 


(3.18) 


(3.19) 


compatible with 2.14[ and differentiation of (|3.6[ ) with respect to x. We use the 
initial data S'(O) = 0 and A1(0) satisfies (2.6) cf. [3H Section 2]. If the matrix 


Af(s) satisfies ( |2.6[ ), then for all s the matrix Al(s) inherits the properties of 

M(0); 


M{s)x{s) = p{s) M{s) = M{sy 


(3.20) 


and 3A4(s) is positive definite on the orthogonal complement of x(s), cf [40] . 
In order to write down such an ip{x) we need to be able to write s as a function 
of the coordinate variables x G M. We know x(s) traces out a smooth curve 
7 in and if we assume x{s) is non-trapping then this curve is a straight 
line when s is sufficiently large. We consider R large enough so that the set 
{x : |x|g/ < 6R} contains M. We set 

fl{eo) = {x : \x\g' < 6R dg>{x,j)<eo} (3.21) 


as a tubular neighborhood of 7 with radius eg in the ball {|x|g' < 6R}. Choosing 
eg sufficiently small we can uniquely define s = s(x) for all x G n(eg) such that 
x(s) is the closest point on 7 to x provided 7 has no self-intersections. The 
variable s is the analogue of the time variable for time dependent problems. 

We now define a cutoff function (p\(x) G C°°(R‘^) for A > 0 such that 


(px{x) 

One can arrange that there 


_ r 0 if xG fl=(A-i/ 2 ^) 
“ \ 1 if xG fI(A-i/ 2 ‘^) 

is a constant C such that 

sup < CA-™. 

xGM 


(3.22) 


(3.23) 


We drop the subscript A for the rest of this paper. We used the fact that 
A“^ < eg so the definition of the parameter s makes sense as above. 

We also recall the following result: 

Lemma 2 (|52 Corollary 5]). Let ip{x) be the phase function of a zeroth order 
beam. We have 


exp{—2XQijj{x)) ~ exp(—AC'dg'(x,x(s))^), 
where C is independent of X. 
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Let B denote the set 


B = {x : dg'{x,x{s)) > \ ^2 0<s< 6 i?}, cr > 0, ct G K. 

We conclude that since 25^/^(s,a:) ~ dg>{x,x{s))‘^, exp(—2A3V'(2^)) is exponen¬ 
tially decreasing in A for all x G B. Notice that we are taking more care to 
construct the cutoff functions than in [52], as they are crucial for the phaseless 
measurements. 

Proof. We need only observe that Af(s) is a bounded and positive definite 
matrix. From the form of the phase functions constructed the desired result 
follows. □ 

The construction of the localized cutoff hnishes the construction of C/a- A 
standard argument gives that U\ extends across each local coordinate chart to 
cover the ray x{s) iteratively, c.f. [551 Section 7]. □ 


4 Introduction of the Source Terms 


We now introduce source functions. We summarize the results from [32] to show 
we can build such sources. We claim: 


Theorem 5. There exists a Gaussian beam solution which solves the equa¬ 
tion (1.2) and is found by solving the Dirichlet problem on one side of hyper¬ 
planes which contain a source point. 


This argument is similar to [52] Section 2.1] and is repeated for completeness. 
We let p be a function such that jVgpj = 1 on the hypersurface S = {cc : p{x) = 
0}. Let Xq be a point in E and we let (a;(s),p(s)) be the solution path — 
in other words the null-bicharacteristics with (a;( 0 ),p( 0 )) = (xg, n(xo)Vgp(xo)). 
The hypersurface S is given by s = a(y) with (t( 0) = 0 and Vct(O) = 0 where x = 
(s,y) and y = (yi,..., yd-i) is transversal. We let the optics Ansatz l/(x) have 
initial data (x(0),p(0)), and be defined in this tubular neighborhood. We let 
(7~^ be the restriction of (7 to {x : p(x) > 0}. Because we need to have a source 
term which is a multiple of 5{p) , we also need a second, ‘outgoing’ solution U~, 
defined on {x : p{x) < 0}. It is then equal to [/+ on the hypersurface E. For 
example, we can write the ingoing and outgoing optics solutions as 


U~^ = A'^{x,X)exp{iX'ijj^{x)), U =A {x, X) exp{iXtp {x)), (4.1) 


where [52] have set = 'ij:~ and A+ = A~ on E. The requirement that their 
Taylor series coincide on the boundary is equivalent to setting dy\y=g'il)^{(j{y)) = 
dy\y^g'tf~{(j{y)). We extend [/+ to be 0 on {a; : p{x) < 0} and U~ to be 0 
on {x : p{x) > 0}. We define our geometric optics Ansatz solution U to be 
U = U~^ -\- U~. We set A = A~^ = A~ on E. In order to add the source term. 
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we notice that 


LU = ix(^ A{x, A) + ^(a;, A) - ^^{x, A)^ exp{iX(j)+)6{p)) 

(4.2) 

+ fgb = 9oHp) + fgb, 

where i'{x) = Vp{x) is the unit normal to E. We consider the singular part 
go5(p) = h(x, X) to be the source term and fgb the error. Then we obtain 

i i 

/ = exp(iAV'^) ^ c+(a::)A“'^ + exp(iAV'“) ^ c~{x)X~\ (4.3) 

i =-2 i =-2 


where we extend to be zero for p{x) < 0 and c^- to be zero for p{x) > 0 . 
We know by construction that c ^2 = ^{dg{x,x{s))^) and = 0{dg{x,x{s))), 
respectively for 1 ®* order beams. 


By construction we see that h{x,X) as defined in (2.2) satisfies (4.2), with 
p{x) a boundary defining function. The incoming and outgoing Gaussian beams 
are matched at the boundary of M by construction. 


Remark 6. While this procedure may seem backwards and a bit ad hoc, it is 
useful for deriving good error estimates which are needed for the use of the 
embedding theorems. 


5 Error Estimates 

We claim for N'*'^ order Gaussian beams: 

Lemma 3. We have the following estimate for the Gaussian beam error: 

(5.1) 

with C independent of X depending on then norm of both the initial 

data and n^(x). 

Proof. We have that the Cj{x)^ are bounded and 

cf{x)= d^^jdg>{x,x{s)f, j =-2,1, (5.2) 

\p\=N-2j-2 

where the dpj are bounded by Taylor’s theorem. We obtain 

|c^(a:)|g' < Cjdgfx,x{s))^~'^^~'^\(t){x)\g^, (5.3) 

with the Cj uniformly bounded independent of A. As Al(s) is a matrix with 
positive definite imaginary part, 

lmf}^{x) > Cdgfx,x{s))^ (5.4) 
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for some positive constant C. From the elementary inequality [32j for a,b > 0, 


If exp(—ab^) < Cptt exp(—ab^/2), 



(5.5) 


with p = N — 2j — 2, a = Xc and b = dgi{x,x{s)), x S Q{X we obtain 


\fgbix)\g> <exp{-XQip^{x)) \cj{x)\g,X f (5.6) 

i =-2 

see M- It follows that 

i 

\fgb{x)\g' < exp{-Xcdg<{x,x{s)f) Y Cjdgi{x,x{s))^~'^^~'^\~^ < 

i=-i 

i 

Cexp{—Xcdgi{x,x{s)Y/2) Y^ <Cexp{—Xcdg'{x,x{s)Y/2)X~^^'^'^^ 

i =-2 

with C a generic constant independent of A depending on then norm 

of both the initial data and n^{x). If we differentiate c^(x) and exp(iXip(x)), 
then we have for the order beams 

ll/sb(2^)lli/™(|x|g,<i?) — (^■'^) 

^^-N+2+2m J exp{- 2 XCdg{x,x{s) f)dx < C;^-W+2+(l-d)/2+2m^ 

as claimed since n(A“^'^^‘^) C n(A“^/^) and the integrand is positive. □ 


Remark 7. Notice that this error now includes a function which is localized in 
a neighborhood ofO{l/V^), which is more precise than in VJiHf . The reason for 


the precise (fx is twofold: We need a localized solution for (8.18) to hold, and 
we are also considering a bounded domain. 


6 Extension of the Ansatz 

The extension of the Gaussian beam Ansatz to M' will provide us with the 
means to make observations further away from the source and still identify 
source terms. In the last section we built an Ansatz which is a subset of a 
slightly extended version of M. We view the manifold M as an embedded 
subset of so we have the following inclusions: 

M cW C {x : \x\g^ <6R}. (6.1) 

for R sufficiently large. We briefly review the results of (35]. We set 

fu = LU-goSip). (6.2) 
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We would like this function to be supported in \x\gi < 6R and be 0(A“^). We 
let G\{x) be the Green’s function for the Helmholtz operator L. In order to 
extend our approximate solution U, we introduce a smooth cutoff Xa(x) such 
that 


Xa{x) = 1 for \x\g> <{a- 1)R, 

Xa{x)=0 for \x\g' > aR. 

Now we set 

U{x) =X3{x)U{x)+ [ Gx{x-y)x 5 {y)L[{i-X 3 {y))U{y)]dy. 


In 


they prove 


U-U 


Hrr^i\x\,<6R) 


= C>(A-”). 


(6.3) 


(6.4) 


Therefore the size of U and the extension depends on the size of fgb- Further¬ 
more, using ideas in Vainberg [50], they also prove 


The triangle inequality allows us to conclude 


u-U 


H^i\x\,<6R) 


< CX ll/gb|li/m(|a;| <i^) 


(6.5) 


( 6 . 6 ) 


The proofs done in [32] are for when (M, g) has the Euclidean metric, but can 
be easily extended as long as M is simple with respect to g. 

In conclusion, we have that: 

Lemma 4. There exists an extension of the Gaussian beam solution of the 
problem (1.2) to the whole space such that 

V-l 


U-U 




Here we know U is supported in |a;| < 6R and we know u is in 


(6.7) 


with 


g' metric so that for R sufficiently large, 
as u—U 

H^(\x\^,<m) 


u-U 




is essentially the same 


7 Integral Transforms Generated by Hamilto¬ 
nian Flows 

We now consider the problem of recovering a function stably from its integrals 
over the integral curves of a Hamiltonian flow. A related nonlinear problem was 
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considered in [5], and generalizations of the geodesic X-ray transform problem 
have been studied before (see e.g. [giiaill]). However, we are not aware of any 
discussion of the present problem in the literature. 

Let (M, g) be a compact Riemannian manifold with boundary and g: M —>■ K 
a smooth function. We let denote a magnetic field which is a closed 2-form. 
We consider the law of motion governed by Newton’s equation; 

V^7 = y(7) - Vg(7) (7.1) 

We have Y : TM —>■ TM is the Lorentz force which is associated to which is 
the map uniquely determined by 




(7.2) 


Vx € M,^,r] G TxM. The curve which satisfies (7.1) is called an MV geodesic 


as per [Q. The MV stands for Maupertius principle. The equation (7.1 1 defines 
a flow (^4 on TM. which we call an MV flow. Whenever = 0 the flow is the 
potential flow and this corresponds to the Hamiltonian flow, corresponding to 
the Hamiltonian function H : T*M —>• K given by H{x,p) = \gij{x)p^p> -I- q{x) 
in local coordinates. For our purposes it is sufficient to examine potential flows, 
but our results could extend to situations with a magnetic term included, in 
which case the Hamiltonian would be H{x,p) = \p + ajg -I- q{x), Q. = da for fl 
exact. 

We refer to the curves generated by the Hamiltonian flow as i7-geodesics. 
For an MV flow, the energy £(x,^) = \\^\^g + q{x) is an integral of motion. The 
law of conservation of energy says for all MV geodesics the energy is constant 
along the geodesic. 

The MV flow depends on the choice of energy level. Consider an i7-geodesic 
7 : [0,T] —>• M. The value of is independent of t, and we write Hq 

for this constant. We assume that Hq > maxg so we can set ShqM = £~^{Ho), 


x^M 


the bundle with energy Hq over M. It is then necessary for Hq to be larger than 
maxg(x). Otherwise we get a.s x G M any vector ^ G ShqM has non-positive 

xGM 

length. We let i/(x) the inward normal to dM at x. We let the set of vectors 


O+ShoM = {(x,e) GShoM: > 0} (7.3) 

For (x,^) G ShqM let t(x,^) be the time such that the MV geodesic with 
7 ( 0 ) = x, 7 ( 0 ) = f reaches dM. The function t{x,^) ■ ^Ho M —M is smooth by 
Lemma A.5 of [5]. 

Let A(-, •) denote the 2"'’* fundamental form of DM. The boundary of dM is 
strictly MV convex if 


A(x, 0 - d^q{v{x)) 

for all (x ,0 e d+SnoM. 

For all X G M, the MV exponential map at x is the partial map 

exp^^ : T^M -G M 


(7.4) 


(7.5) 
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which is given by the following: 

= IT o t>0,{x,^)£SHoM (7.6) 

For all X G M, exp^^ is a smooth partial map and T^M is C°° smooth on 
TxM/{0}. We let exp^ denote the exponential map associated to the potential 
Hamiltonian (H = 0). In analogy to the usual exponential map exp^: T^M —> 
M, we have the Hamiltonian exponential map exp^: T^M —>■ M for each x G 
M, defined by the Hamiltonian flow. 

Since our manifold has a boundary, this is only a partial map. We say that 
the manifold is simple at energy Hq if the map exp,^: {exp^)~^{M) -G M is a, 
diffeomorphism for each x G M, and also that dM is strictly j\4P convex. We 
remark that simplicity implies ffg > maxq. It is more stringent to assume a 

xGM 

system to be simple than to assume it to be non-trapping. 

We define the Hamiltonian flow transform of a function /: M —>■ K at energy 
iJo for (a;,^) G B+ShoM as 


iHofh) = f fh(t))dt, (7.7) 

^0 

IhJ :C°°iM)^C°°{d+SHoM) (7.8) 

where 7 is an i7-geodesic of energy Hq. We denote the corresponding normal 
operator by with the adjoint operator 

I*hJ : C^{d+SHoM) -G C~(M). (7.9) 

For a fixed energy level Hq we let 7 (t) be an AiV geodesic. We consider the 
change of variables 


t 

s{t)= j 2{HQ-q{^{t))dt. (7.10) 

0 

It follows that s is the arc-length of cr(s) = 7(t(s)) under the metric g = 2{Hq — 
q)g. A version of Maupertuis principle says that cr(s) = 7(^(5)) is a unit speed 
magnetic geodesic of the system ( 5 , 0 ). We quote these results from [B] without 
proof. 

Theorem 8 (Theorem 2.1 in [^). Let {M,g,q) denote an AiV potential system 
on M and Hq a constant such that Hq > maxg(a;). We suppose ^{t) is an 

x£M 

MP geodesic of Hq, then a{s) = 7(t(s)) is a unit speed geodesic of the system 
(M, g, 0) on M. 

and 

Proposition 1 (Prop 2.2 in [6]). The AiV potential system {M,g,q) on M of 
energy Hq is simple if and only if so is system {M, g, 0). 
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We use these two results to prove the following: 

Lemma 5. Suppose (M,g,q) is simple at energy Hq w.r.t. the Hamiltonian 
H. Define a new metric g conformal to g by g = 2 {Hq — q)g. Then {M,g) 
is a simple Riemannian manifold and a reparametrization turns H-geodesics of 
energy Hq to unit speed geodesics w.r.t. g. 

Moreover, if I denotes the X-ray transform w.r.t. the metric g, then Inofil) = 
I[f /2 {Hq — 9)](7 o r), where r is the reparametrization. 

Proof. The first part follows from 0 Proposition 2.1] and the second part is a 
straightforward calculation. □ 

Let us remark that the boundary dM is strictly convex with respect to g if 
and only if A(u, v) > —di,q{x) for all {x, v) G with x G dM and v _L dM. 

Here di, is the inward normal derivative and A(-, •) the second fundamental form. 
This calculation can be found in [HI Lemma A.4]. 

Theorem 9. Suppose dimiW > 2. If the manifold is simple for Hq, then Ihq 
is injective on L‘^[M). Moreover, we have the stability estimate 

where M' G M is a slightly extended manifold and C is a constant depending 
on the manifold, q and Hq . 

Proof. By Lemma[^it suffices to show that the X-ray transform on the manifold 
{M,g) is injective and stable. Injectivity [131 Theorem 7.1] and stability |43l 
Theorem 3] are known for simple manifolds. □ 

In fact, it suffices to assume that the underlying manifold {M,g) is simple if 
Hq can be taken arbitrarily large. The following corollary follows immediately 
from Theorem 1^ above and Lemma below. 

Corollary 2. Let {M, g) be a simple manifold of dimension two or higher. Then 
for sufficiently large Hq 


II/IIl2(m) ^ C\\NHof\\HPM-)^ (7-12) 

where M' G M is a slightly extended manifold and C is a constant depending 
on the manifold, q and Hq. 

The corollary applies, in particular, to the closure of any strictly convex and 
bounded Euclidean domain. 

Lemma 6 . If {M, g) is simple, then the manifold M is simple at energy Hq for 
large enough Hq. 

Proof. It is well known that simplicity is an open condition: small perturbations 
of simple metrics are still simple. Therefore for Hq large enough the metric 
(1 — q/HQ)g is simple. Rescaling the metric does not alter simplicity, so also the 
metric g = 2 {Hq — q)g is simple. The claim then follows from Lemma□ 
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Remark 10. This remark may be of interest for integral geometers, although 
not strictly relevant for this paper. 

In a similar way we can prove injectivity results for the Hamiltonian ray 
transform of tensor fields. For a smooth tensor field f of order m the formula 
of lemma [5^ is iHofil) = — 2q)"^ ^/](7 o ’’)■ If solenoidal injectivity is 

known for the geodesic X-ray transform of tensor fields on the manifold {M,g), 
then we know that (2 Hq — f = adh for some tensor h of order to — 1 

(here dg is the covariant derivative w.r.t. the metric g and a is symmetrization). 
On simple surfaces solenoidal injectivity is known for tensors of all orders \Ji(A 
Theorem 1.1], and there are several solenoidal injectivity results for higher di¬ 
mensions as well (see eg. In particular, solenoidal injectivity is 

known for to = 0 and m = 1 (see W)- We refer to the survey article for 
more details on known results. 

Ifm= 1, I Hof — 0 and solenoidal injectivity of I imply that f = dh for some 
smooth scalar function h. If m > 2, the gauge condition is more complicated. 
We will not pursue this direction further here, as it is irrelevant for our main 
problem. 


8 Observability Inequalities 


We prove the main theorem, Theorem[2in this section. We consider our globally 
defined complex optics solutions Ui and C/ 2 , which were constructed previously. 
Dropping the superscripts xq, wq where it is understood, from our approximation 
we know that the main term of interest is 


|Ml| 


F 2 II = 


\U 1 \-\U 2 \ +o(a-i/2) 


( 8 . 1 ) 


for X € dM. This estimate is a result of building an approximate solution with 


N sufficiently large. We conclude from (6.6) and (5.7) 


Ui - ui 


< 


C 


H-iM') y/X 


( 8 . 2 ) 


Here C is a generic constant, which depends on the norm of nf(x) 

with N such that max{(l — d)/2 + 2 + 2s, 0} + 1. We use the fact that Ui, and 
ui are bounded in norm to conclude 


Ui - Ml 


< 


c 


CO(M') yX 


(8.3) 


by Sobolev embedding, provided s > d/2 

We use the estimate on the zeroth order terms 


\Ui — exp(iA'i/’) 


I ^ 

lco(M') - y- 


(8.4) 


where is the first order term, found to leading order by (3.2) and (3.18). 
Combining the estimates (8.4) and (8.3), gives the estimate (8.1). Now we need 
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to combine the estimates to recover the X-ray transform. We would like to use 


the first order coefficients to reconstruct the coefficient n^{x). Using (3.2), we 
know that 


a^{x) exp{iXipi{x)) — 09 ( 0 ) exp I — / {n1{xi{t)) dt — J trAf i(f) dt 1 exp{zA'!/'i(a;))^(x) 

Vo 0 / 

(8.5) 


by choice of the cutoff function in (3.22). 


Examining (8.1), we are interested in the left hand side. Using (8.5), we 


would like to approximate it by 


ao(0) ^exp ^ J nl{xi{t))dtj - exp ^ J ril{x 2 {t)) ^ 

exp trA4i{t)dt \ exp{\tpi{x))(l){x) , 


( 8 . 6 ) 


where 


'01 (x) = —(x — xi(s))QA4i(s) ■ (x — xi(s)). 


(8.7) 


We need the following lemma to control the error, which works for level sets 
of general Hamiltonians H, however we chose to focus on the case H{x,p) = 
\p\g ~ which is relevant for our solution curve. Again, we remind the reader 
the equation ( 3.17[ ) is the correct one that we want to solve to high order by the 
derivation of Lemma Recall Grownwall’s identity: 

Theorem 11 (Grownwall’s Inequality). Let w : [ffiAi] be continuous 

and non-negative and suppose that w obeys the integral inequality 


I 

'}{t) < A-\- J B{s)w{s) ds 


for all t G [toUi] whenever A > 0 and B : [toUi] 
non-negative. Then we have 


']{t) < A exp I J B{s) I ds 


( 8 . 8 ) 

is continuous and 

(8.9) 


for all t G [to) ti]- 

We let e G (0,1) be a fixed positive constant. We set H{x,p) = p-\/pH{x,p). 
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Lemma 7. (Uniqueness) Assume that ipi,ip 2 € C°°{M;C), sq^s G {a,b), satisfy 

fji{x) = il)2{x) + e>(((ig'(a;,a:i(so))”"+^), for sq {da;'tpj){xj{s)) =Pj{s), 
dsipjix) + H{x,Wgi)j{x)) = 

for some m € N. Then the following holds 

sup |a;i(s) - a; 2 (s)| = ^(e) sup \Mi{s) - M 2 {s)\ = 0{e) (8.10) 

s^{a,b) s^{a,b) 

Proof. We can map the neighborhood of the flow on the manifold to the Eu¬ 
clidean plane and use the metric gij = Sij the standard Euclidean one. In this 
case we have by definition of the Hamiltonian flow. 

1 ^ < < \ini -n 2 )(a;i(s))| -bClxl (8.11) 

where C = sup(Vn 2 (a;)) by the mean value theorem and x(s) = a;i(s) — a; 2 (s). 

X 

Given initial conditions x(0) =0 we can use GronwalFs inequality with w = \x\ 

t 

and A = f eels to conclude 
*0 


\x{s)\=0{e) ( 8 . 12 ) 

since | |nf — ^il Ic 3 (m') ^ ^ terms depend on the length of the geodesic. 

We also have by definition of the matrix Ricatti equations: 

< l^■^(s)l = -D^nlixiis)) + D^nlix2is)) + MUs) - MI{s) 

(8.13) 

which then implies 

|-^Ai(s)| < - n 2 )(a;i(s))| -f Ci\x\ -b C' 2 |M| (8.14) 

ds 

where Ci = sup(Zl^n 2 (a;)) by the mean value theorem and M{s) = A4i(s) — 

X 

A42{s), C 2 = sup {A4j(s)}. Applying Gronwall’s theorem again gives the 

1=1, 2 ,s 

second desired inequality. □ 


The expansion implies 

(x-a;i(s)) = (a;- 0 : 2 ( 5 )) +0{e) + 0{dg>{x,xi{s)f); (8.15) 

(x — xi(s))9AIi(s)(x — xi(s)) = 

(x — X2(s))9Al2(s)(x — ^ 2 ( 5 )) -b 0(e) -b 0{dg'(x, xi{s))^) 

and also 


AMi(s) = 9A42 (s) -bO(e) 


(8.16) 


20 



It follows Lemma [7] in local coordinates: 


\U 2 {x)\ =|ao(0)exp y J dtj exp y J exp(AV’i(a:))(/) 2 (a;)|x 

(8.17) 

(exp(-C>(Ae)) 

with V'i(3^) 3.S above. Since Lemma [t] is local, the order terms depend on 
diam(M) and C^{M) norm of n\{x) and n\{x), made precise by (8.15),(8.16). 

This allows us to conclude (8.6), with appropriate loss of error. Notice that 


even if e > if we can re-define U 2 with a 4 ‘ 2 {x) which has larger significantly 
larger support and still obtain the same estimate as (8. since by Lemma 1^ 
the contribution of the error estimates outside of a neighborhood of 
is negligible. The cutoffs were kept with support in to illustrate that 

the integration in a neighborhood around the tube is the most important con¬ 
tribution, c.f. m- 


We now examine (8.6). Because 5Ali(s) is a positive definite matrix, by 
Corollary we have 


sup 

xeM 


eyip{XilJi{x))(j){x) 


= 1 . 


(8.18) 


Indeed, xi(s) reaches the boundary and is contained in Q{X ^/^). Because the 


other coefficients of (8.6) are independent of x, we obtain that the supremum 


of (|8.6|) over x S dM equals 

t(xo,ujo) 


t{xo,uio) 


C{s) 


Here 


a()(0) I exp I — J a{n\{xi{t)) dt j — exp | — J Q;n 2 (x 2 (t)) J x E{x) 

(8.19) 


/ r(a:o,‘^o) 


C{s) = exp 


trAli(t) dt 


( 8 . 20 ) 


V 


with t(xo,ujo) the time it takes for x(s) defined by (2.5) with initial conditions 
(xo,cuo) to exit M. We also have 


\E{x)\\co(m') = exp{0{-Xe)) 


( 8 . 21 ) 


by (8.17) from Lemma (0. 

We need the following lemma to obtain the ray transform which is taken 
from [H]. 
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Lemma 8. Let A(x) and B{x) he positive functions in and e G (0,1) 

such that 


||exp(-A(a;)) - exp(-B(a;))||^o(R) < £■ (8.22) 

Then there is a constant C depending on the (^^(K) norms of A and B, such 
that 


P(a:) -B(x)||co(r) < Ce. 


(8.23) 


Proof. By the mean value theorem, there exists an r* between B{x) and A(x) 
for each fixed x such that 


|(exp(-P - exp{-B))\ = 


' f \ 

— / exp(—r) dr 


= |i? - ^1 exp(-r*). (8.24) 


The desired result follows by taking the supremum over x and applying (8.22). 


□ 


Now we can see by assumption that 5, A ^ are less than eq- Applying Lemma 
([^ with e = eg and = A, Iho '^2 = B we obtain: 


’^ 2)1 lc’0(a5M+) — \/S' ) • 




(8.25) 


Cl denotes a generic constant depending on the C^{M) norm of nj, n|, while 
C 2 depends on the norm. Hq denotes the flow associated to n^x). 

Here we have used the fact that we can Taylor expand ( 3 : 2 ( 5 )) around a:i(s), 
and a:i(s) — 0 : 2 ( 5 ) = 0(e) via ( |8.15 ) and the smallness condition (2.8). 

Now we set n\ — n\ = h?, and we recall that it has compact support. Because 
we have assumed that M is strictly convex, it follows from (8.25) that 


|^ffo(^l 1(71(95^+) — 


C2 

A/3' 


(8.26) 


by the embedding theorems in the Appendix with fc = 2, t = 1. Here Ci denotes 
a generic constant depending on the C‘^{M) norm of nf, n|. The compact 
support of hf implies that 


\ri 


Ha(Pi ^ 2 ) I lci(95M+) —^1 i ^ 


Ch 

A/3' 


(8.27) 


To finish the proof of Theorem[^ we use (8.27) along with the stability estimate 
of corollary!^ In particular, there is a number c' > 0 that depends on M, g 
and ni such that the following holds. If we choose initial speeds with |a;o| = c', 
then the flow transform is stably invertible. We notice that nf G C^{M) and 
also by choice of nf = = 1 on dM the stability estimate over M' makes no 

difference, so that the Corollary [^applies with n^{x) = —q{x) and Hq = 0. 
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9 Appendix: Embedding theorems 

This result is taken from Section 2.4 in [^. We let /i, /2 G C^{M) with fc > 1. 
We assume there is a constant C > 0 and an e << 1 such that 


ll/l|lc'=(M) ”1” II/2 |Ic'=(M) — ^ ll/l f2\\c(M) — ^ (9-1) 


We use standard interpolation estimates, see Section 4.3.1 in |49] to normalize 
the use of norm estimates. As an example using interpolation estimates we 
obtain for generic / G C^{M) with fc > 1 



(9.2) 


with 0 < 0 < 1 with ti > 0,^2 > 0. As an immediate consequence of the 
assumptions on fi , /2 we have that 


\\fl-f2\\cHM)<C^'^ 


(9.3) 


for each t > 0 if fc > t. 
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